Getting started with geospatial data

• What are we doing today?
  • Shapefiles
    • Polygons
    • Points
    • Lines
    • Mapping with the package sf
  • Coordinate reference systems
    • Latitude/longitude
    • Projections

Shapefiles

• Shapefiles are a common format for geospatial data
  • They are a form of vector data
• Shapefiles are made up of at least three files:
  • .shp - the shape itself
  • .shx - the index
  • .dbf - the attributes
  • .prj - the projection
    • This one is not technically necessary! But it’s common to have.
  • What these all mean isn’t important for now, just make sure they are there! Check the day2files folder on github.

Let’s look at Northern Malawi

Types of features

• Polygons
  • Areas
  • Districts, countries, etc.
• Lines
  • Lines
  • Roads, rivers, etc.
• Points
  • Points

Let’s start with polygons

Imagine a rectangle, on a coordinate plane

Imagine a rectangle, on a coordinate plane

Imagine a rectangle, on a coordinate plane

All features are made of vertices

All features are made of vertices

One more example of polygons

One more example of polygons

Lines

• Lines are also made up of vertices
• But they are not closed

Lines example - “Primary” roads in India (2014)

One road

Length of this line feature: 1619.93 (m)

Points

• Points are exactly what they sound like: points!

• What could be a point?

Points

• Points are exactly what they sound like: points!

• What could be a point?

  • A city
  • A weather station
  • A tree
  • A household
  • etc.

What do you think this is?

What do you think this is?

Train stations! First 15

The train stations are a collection of features

• Just like before, the shapefile is a collection of features!
• The only difference now is that each feature is a point

Reading shapefiles in R

• My go-to package for shapefiles in R is sf
• Reading shapefiles is VERY easy! And you can treat them like dataframes.

library(sf)
# this is the shapefile for Gambia
gambia <- read_sf("day2files/gambia.shp")
gambia

Simple feature collection with 115 features and 17 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: -16.81426 ymin: 13.04347 xmax: -13.79112 ymax: 13.82274
Geodetic CRS:  WGS 84
# A tibble: 115 × 18
   Shape_Leng Shape_Area ADM3_EN       ADM3_PCODE ADM3_REF ADM3ALT1EN ADM3ALT2EN
        <dbl>      <dbl> <chr>         <chr>      <chr>    <chr>      <chr>     
 1     0.0115 0.00000688 New Town East GM010101   <NA>     <NA>       <NA>      
 2     0.0206 0.0000161  New Town West GM010102   <NA>     <NA>       <NA>      
 3     0.0200 0.0000192  Soldier Town  GM010103   <NA>     <NA>       <NA>      
 4     0.0337 0.0000302  Box Bar       GM010201   <NA>     <NA>       <NA>      
 5     0.156  0.000591   Campama       GM010202   <NA>     <NA>       <NA>      
 6     0.0179 0.0000127  Crab Island   GM010203   <NA>     <NA>       <NA>      
 7     0.0322 0.0000533  Half Die      GM010301   <NA>     <NA>       <NA>      
 8     0.0216 0.0000236  Jollof Town   GM010302   <NA>     <NA>       <NA>      
 9     0.0203 0.0000257  Portugiese T… GM010303   <NA>     <NA>       <NA>      
10     0.494  0.00937    Basse         GM020101   <NA>     <NA>       <NA>      
# ℹ 105 more rows
# ℹ 11 more variables: ADM2_EN <chr>, ADM2_PCODE <chr>, ADM1_EN <chr>,
#   LGA_EN <chr>, ADM1_PCODE <chr>, ADM0_EN <chr>, ADM0_PCODE <chr>,
#   date <date>, validOn <date>, validTo <date>, geometry <POLYGON [°]>

Plotting is also very easy

ggplot() + 
  geom_sf(data = gambia)

My go-to theme

ggplot() +
  geom_sf(data = gambia) +
  theme_bw()

Other changes you can make

ggplot() +
  geom_sf(data = gambia, fill = "#f0f1eb", color = "#006334") +
  theme_bw() +
  labs(subtitle = "Admin3s in Gambia")

Give it a try with TAs (admin3) in Malawi (mw3.shp)

library(sf)
# this is the shapefile for the northern region of Malawi, TA level
northmw <- read_sf("day2files/mw3.shp")
ggplot() +
  geom_sf(data = northmw)

ggplot() +
  geom_sf(data = northmw) +
  theme_bw() +
  labs(subtitle = "TAs in Northern Malawi")

One more example - map from earlier

admin2 <- read_sf("day2files/mw2.shp")

ggplot() + 
  geom_sf(data = admin2, 
    fill = "white", color = kdisgray) +
  geom_sf(data = admin2 |> filter(DIST_CODE=="107"), 
    fill = "white", color = kdisgreen) +
  theme_bw() +
  labs(subtitle = "Districts in Northern Malawi")

What if we want to plot mw2, mw3, and mw4?

What if we want to plot mw2, mw3, and mw4?

admin2 <- read_sf("day2files/mw2.shp")
admin3 <- read_sf("day2files/mw3.shp")
admin4 <- read_sf("day2files/mw4.shp")

g1 <- ggplot() + 
  geom_sf(data = admin2, fill = "white", color = "black") +
  theme_bw() + labs(subtitle = "A. Admin2 (districts)")
g2 <- ggplot() + 
  geom_sf(data = admin3, fill = "white", color = "black") +
  theme_bw() + labs(subtitle = "B. Admin3 (TAs)")
g3 <- ggplot() + 
  geom_sf(data = admin4, fill = "white", color = "black") +
  theme_bw() + labs(subtitle = "C. Admin4 (EAs)")

What if we want to plot mw2, mw3, and mw4?

• Enter cowplot!

library(cowplot)

plot_grid(g1, g2, g3, ncol = 3)

A quick note about shapefile sizes

• What do you think is the main determinant of the size of a shapefile?

And that’s only shapefile

• Other geospatial data can get even bigger!

• How large do you think the folder on my computer that contains imagery for all of Malawi (at 5m resolution) is?

Latitude and longitude on a globe

• The most common coordinate reference system (CRS) is latitude/longitude
  • Latitude: North/South
  • Longitude: East/West
  • The equator is at 0° latitude
  • The prime meridian is at 0° longitude
    
    
• But there’s a problem with using latitude/longitude
  • The Earth is a sphere (well, more or less; really an oblate spheroid)

The basic problem

• The basic problem is that one degree of longitude changes at different latitudes!
  • At the equator, one degree of longitude is about 111 km
  • At 15N/S, one degree of longitude is about 107 km
  • At 30N/S, one degree of longitude is about 96 km
  • At 45N/S, one degree of longitude is about 79 km
  • At 60N/S, one degree of longitude is about 56 km
    • This explains Greenland!
      
      
• It’s not an easy problem to solve, as all solutions have drawbacks!

Preserve shape, give up area

Projections

• A projection is a way to represent the Earth’s surface on a flat plane
  • In other words, it’s a way to transform the three-dimensional surface of the earth into two dimensions
  • Think of this as a way to “flatten” the Earth

Robinson projection: compromise across the board

Meractor projection: preserves angular relationships

Azimuthal Equidistant projection

Mollweide Equal Area Cylindrical projection

Three families of projections

• The QGIS website has a great explanation of the three families of projections
  • Cylindrical
  • Conical
  • Planar

Three families of projections

Coordinate reference systems

• Coordinate reference systems (CRS) are a way to define how coordinates are represented
  • This includes the projection, but also other things
• The most popular is probably WGS 84 (EPSG:4326), which is a geographic CRS
  • This is latitude/longitude
• One degree of lat/lon:
  • 60 minutes
• One minute:
  • 60 seconds

Geographic vs. projected

A common CRS: Universal Transverse Mercator (UTM)

• Universal Transverse Mercator (UTM) is a projected CRS
  • It divides the world into 60 zones
  • Each zone is 6° wide
  • The equator is the origin of each zone
  • The equator is at 0 m

Universal Transverse Mercator (UTM)

A common CRS: Universal Transverse Mercator (UTM)

• UTM defines X values (“longitude”) FROM THE MIDDLE of each zone

  • This middle line is called the central meridian

• UTM defines Y values (“latitude”) from the equater

• There are some other details:

  • UTM values are never negative, so we offset values
    • This is called a “false easting” and “false northing”
    • Details aren’t super important

• Note: it ignores altitude. It assumes a perfect ellipsoid

“Transverse mercator” projection per zone (N and S)

Projections with sf

admin2 <- read_sf("day2files/mw2.shp")

# Geographic CRS: WGS 84 - lat/lon
crs(admin2)

[1] "GEOGCRS[\"WGS 84\",\n    DATUM[\"World Geodetic System 1984\",\n        ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n            LENGTHUNIT[\"metre\",1]]],\n    PRIMEM[\"Greenwich\",0,\n        ANGLEUNIT[\"degree\",0.0174532925199433]],\n    CS[ellipsoidal,2],\n        AXIS[\"geodetic latitude (Lat)\",north,\n            ORDER[1],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        AXIS[\"geodetic longitude (Lon)\",east,\n            ORDER[2],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n    ID[\"EPSG\",4326]]"

What is the appropriate zone for Malawi?

• Go to Google
  • Search “UTM CRS Malawi”
• What did you find?

What is the appropriate zone for Malawi?

admin2proj <- st_transform(admin2, 20936)

# It has changed! Now it's projected
crs(admin2proj)

[1] "PROJCRS[\"Arc 1950 / UTM zone 36S\",\n    BASEGEOGCRS[\"Arc 1950\",\n        DATUM[\"Arc 1950\",\n            ELLIPSOID[\"Clarke 1880 (Arc)\",6378249.145,293.4663077,\n                LENGTHUNIT[\"metre\",1]]],\n        PRIMEM[\"Greenwich\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        ID[\"EPSG\",4209]],\n    CONVERSION[\"UTM zone 36S\",\n        METHOD[\"Transverse Mercator\",\n            ID[\"EPSG\",9807]],\n        PARAMETER[\"Latitude of natural origin\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8801]],\n        PARAMETER[\"Longitude of natural origin\",33,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8802]],\n        PARAMETER[\"Scale factor at natural origin\",0.9996,\n            SCALEUNIT[\"unity\",1],\n            ID[\"EPSG\",8805]],\n        PARAMETER[\"False easting\",500000,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8806]],\n        PARAMETER[\"False northing\",10000000,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8807]]],\n    CS[Cartesian,2],\n        AXIS[\"(E)\",east,\n            ORDER[1],\n            LENGTHUNIT[\"metre\",1]],\n        AXIS[\"(N)\",north,\n            ORDER[2],\n            LENGTHUNIT[\"metre\",1]],\n    USAGE[\n        SCOPE[\"Engineering survey, topographic mapping.\"],\n        AREA[\"Malawi. Zambia and Zimbabwe - east of 30°E.\"],\n        BBOX[-22.42,30,-8.19,35.93]],\n    ID[\"EPSG\",20936]]"

They look quite similar!

They look quite similar! Different coordinates

g1 <- ggplot(admin2) +
  geom_sf(fill = "white", color = kdisgreen) +
  theme_bw() +
  labs(subtitle = "A. Geographic CRS") +
  theme(plot.background = element_rect(fill = "#f0f1eb", color = "#f0f1eb"))
g2 <- ggplot(admin2proj) +
  geom_sf(fill = "white", color = kdisgreen) +
  theme_bw() +
  labs(subtitle = "B. Projected CRS") +
  theme(plot.background = element_rect(fill = "#f0f1eb", color = "#f0f1eb")) +
  coord_sf(crs = st_crs(20936))

Malawi is small. What about something larger?

Now it’s your turn!

• Use the Indian roads shapefile (indiaprimaryroads.shp)

• Do the following:

  • Find the CRS of the shapefile
  • Transform the shapefile to UTM zone 44N (EPSG: 24344)
  • Graph both side by side using cowplot

The solution

roads <- read_sf("day2files/indiaprimaryroads.shp")
crs(roads)

[1] "GEOGCRS[\"WGS 84\",\n    DATUM[\"World Geodetic System 1984\",\n        ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n            LENGTHUNIT[\"metre\",1]]],\n    PRIMEM[\"Greenwich\",0,\n        ANGLEUNIT[\"degree\",0.0174532925199433]],\n    CS[ellipsoidal,2],\n        AXIS[\"geodetic latitude (Lat)\",north,\n            ORDER[1],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        AXIS[\"geodetic longitude (Lon)\",east,\n            ORDER[2],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n    ID[\"EPSG\",4326]]"

roadsproj <- st_transform(roads, 24344)
crs(roadsproj)

[1] "PROJCRS[\"Kalianpur 1975 / UTM zone 44N\",\n    BASEGEOGCRS[\"Kalianpur 1975\",\n        DATUM[\"Kalianpur 1975\",\n            ELLIPSOID[\"Everest 1830 (1975 Definition)\",6377299.151,300.8017255,\n                LENGTHUNIT[\"metre\",1]]],\n        PRIMEM[\"Greenwich\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        ID[\"EPSG\",4146]],\n    CONVERSION[\"UTM zone 44N\",\n        METHOD[\"Transverse Mercator\",\n            ID[\"EPSG\",9807]],\n        PARAMETER[\"Latitude of natural origin\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8801]],\n        PARAMETER[\"Longitude of natural origin\",81,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8802]],\n        PARAMETER[\"Scale factor at natural origin\",0.9996,\n            SCALEUNIT[\"unity\",1],\n            ID[\"EPSG\",8805]],\n        PARAMETER[\"False easting\",500000,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8806]],\n        PARAMETER[\"False northing\",0,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8807]]],\n    CS[Cartesian,2],\n        AXIS[\"(E)\",east,\n            ORDER[1],\n            LENGTHUNIT[\"metre\",1]],\n        AXIS[\"(N)\",north,\n            ORDER[2],\n            LENGTHUNIT[\"metre\",1]],\n    USAGE[\n        SCOPE[\"Engineering survey, topographic mapping.\"],\n        AREA[\"India - onshore between 78°E and 84°E.\"],\n        BBOX[8.29,78,35.5,84.01]],\n    ID[\"EPSG\",24344]]"

The solution

g1 <- ggplot(roads) +
  geom_sf(fill = "white", color = kdisgreen) +
  theme_bw() +
  labs(subtitle = "A. Geographic CRS")
g2 <- ggplot(roadsproj) +
  geom_sf(fill = "white", color = kdisgreen) +
  theme_bw() +
  labs(subtitle = "B. Projected CRS") +
  coord_sf(crs = st_crs(24344))

plot_grid(g1, g2)

The solution

More practice

• I’d like you to find a shapefile for a country of your choice
  • Ths isn’t always easy
  • I always use Google to find one
  • One place you can often find something on humdata.org (almost always shows up in the search!)
• Do the following:
  • Find the CRS of the shapefile
  • Transform the shapefile to the appropriate UTM zone for the country
    • Note that many countries have multiple UTM zones! You can just choose one
  • Plot them side by side using cowplot

Returning to the files

• Shapefiles are made up of at least four files:
  • .shp - the shape itself
    • The geometry
  • .shx - the index
    • The “index”. You can actual recover the .shx file from the .shp file
  • .dbf - the attributes
    • The attributes for each features. Could be names, population values, etc.
  • .prj - the projection
    • The projection information, which we just discussed
• There are often more files, but these are the main ones
  • You must have the first three. The fourth is optional, but common.